"""
Collection of algorithms to find an elliptic curve over Q(sqrt5) from a Hilbert
modular form over Q(sqrt5).

Right now this is starting as an implementation of the method proposed by Dembele
in the case that the conductor is not a perfect square.
"""

import math

from sage.all import NumberField, SR, CC, RR, ZZ, prime_range, parallel, fast_callable, pi, sqrt, real, imag,  EllipticCurve_from_c4c6, j_invariant_qexp, Integer, cartesian_product_iterator, gcd, lcm, continued_fraction, IntegerVectors, factor, prod, Set, srange, IntList, PowerSeriesRing, ComplexField

import sage.rings.number_field
from psage.ellcurve.lseries.aplist_sqrt5 import aplist as compute_aplist
from psage.number_fields.sqrt5 import primes_of_bounded_norm
from sage.schemes.elliptic_curves.period_lattice import reduce_tau
from psage.ellcurve.lseries.lseries_nf import lseries_dokchitser, anlist

from ellcurve_pyx import search_possibilities

import time
import sys

x = SR.var('x')                                 # We are only using this to construct our number
                                                # field. Is there a better way?
K = NumberField(x**2 - x - 1, name='a')
OK = K.ring_of_integers()
a = K.gen()

phi2, phi = K.embeddings(RR) # we choose phi to be the embedding that takes a to positive real number
a1 = phi(a)
a2 = phi2(a)

# This is a workaround for an annoying bug...
def custom_ideal_hash(self):
    """
    hashing of number field ideals is currectly broken, so
    we put a stupid workaround here.
    """

    return str(self.pari_hnf()).__hash__()

sage.rings.number_field.number_field_ideal.NumberFieldFractionalIdeal.__hash__ = custom_ideal_hash


# TODO: make this class better. Right now building characters is rather slow, and evaluating them
# doesn't seem to be so great either.
class QuadraticCharacter:
    """
    Quadratic Dirichlet characters for Q(sqrt5)

    This implementation probably isn't very great right now.

    EXAMPLES::

        sage: from psage.modform.hilbert.sqrt5.ellcurve import K, a, QuadraticCharacter
        sage: chi = QuadraticCharacter(3)
        sage: chi(a)
        -1
        sage: chi(1-a)
        -1
        sage: chi(a) * chi(1-a) == chi(-1)
        True
        sage: chi = QuadraticCharacter(3*a - 1)
        sage: chi(2)
        -1
        sage: chi(3)
        1
        sage: chi(6)
        -1
        sage: chi(11)
        0
    """

    def __init__(self, P):
        """
        INPUT:
            - P -- a nonzero prime or prime ideal of the ring of integers of Q(sqrt(5))

        EXAMPLES::
            sage: from psage.modform.hilbert.sqrt5.ellcurve import K, a, QuadraticCharacter
            sage: chi = QuadraticCharacter(3*a - 1)
            sage: chi(2)
            -1
            sage: chi(3)
            1
            sage: chi(6)
            -1
            sage: chi(11)
            0
        """
        P = K.ideal(P)
        self.k = P.residue_field()
        self.K = self.k.lift(self.k(1)).parent().number_field()
        self.character_table = {}
        self.character_table[self.k(0)] = 0
        self._conductor = P

        g = self.k.multiplicative_generator()
        a = g
        current_sign = -1
        for n in range(self.k.order() - 1):
            self.character_table[a] = current_sign
            current_sign = -current_sign
            a = a * g
        # a should now be gen. We use this to make sure we haven't made a big mistake:
        if a != g:
            print a, g
            raise ValueError("something went wrong.")
        

    def __call__(self, a):
        """
        In the case that a is an element of O_K, return this character evaluated at a.

        If a is an ideal, for convenience we return the character evaluated at a totally
        positive generator of a. The character is not actually well defined on ideals, but
        it takes the same value on every total positive generator (since all such generators
        differ by an ever power of a unit), and when we twist an L-function, this is
        what we will, want to twist by.

        EXAMPLES::
            sage: from psage.modform.hilbert.sqrt5.ellcurve import K, a, QuadraticCharacter
            sage: chi = QuadraticCharacter(K.primes_above(11)[0])
            sage: chi(3)
            1
            sage: chi(a)
            -1
            sage: chi(11)
            0
            sage: chi(K.primes_above(19)[0])
            -1
        """
        if isinstance(a, sage.rings.number_field.number_field_ideal.NumberFieldFractionalIdeal):
            a = make_tot_pos(a.gens_reduced()[0])
        a = self.k(a)
        return self.character_table[a]

    def gauss_sum(self):
        """
        Return the gauss sum for this character. (Don't really know right now if this is correct.)
        """
        I = CC(0, 1)
        characteristic = self.k.characteristic()
        x = CC(0)
        
        d = K(K.disc()).sqrt()
        #d = sqrt(K.disc())
        m = self.conductor().gens_reduced()[0]

        for a, value in self.character_table.iteritems():
            z = value * CC(2 * pi * I * (self.k.lift(a)/(d * m)).trace()).exp()
            #print a, value, z
            x = x + value * CC(2 * pi * I * (self.k.lift(a)/(d * m)).trace()).exp()
            #x = x + value * CC(2 * pi * I * self.k.lift(a).trace()/(d * characteristic)).exp()

        return x

    def conductor(self):
        return self._conductor

def Lone_twist(ap_list, sign, level, twist=None):
    """
    Compute the L-function with the given aplist and given sign, twisted by chi, at 1.
    This should be an L-function for an elliptic curve over Q(sqrt5), and we are going
    to assume that the twist has the same sign as the original L-function.

    ap_list should be a list of pairs (prime ideal, ap) sorted in the same
    order as primes_of_bounded_norm()
    """
    
    from psage.number_fields.sqrt5.prime import primes_of_bounded_norm, Prime
    from psage.ellcurve.lseries.helper import extend_multiplicatively
    from sage.libs.lcalc.lcalc_Lfunction import Lfunction_D


    bound = ap_list[-1][0].norm()
    an_list = IntList(ZZ(bound) + 1)
    an_list[1] = 1
    

    # we are just copying some of William's code here...

    bad_primes = set([Prime(P) for (P,e) in level.factor()])

    # We compute the local factors of the L-series as power series in ZZ[T].
    P = PowerSeriesRing(ZZ, 'T')
    T = P.gen()
    # Table of powers of T, so we don't have to compute T^4 (say) thousands of times.
    Tp = [T**i for i in range(5)]

    # For each prime, we write down the local factor.
    L_P = []
    
    if twist is None:
        for P, a_p in ap_list:
            inertial_deg = 2 if P.is_inert() else 1
            if P in bad_primes:
                # bad reduction
                f = 1 - a_p*Tp[inertial_deg]
            else:
                # good reduction
                q = P.norm()
                f = 1 - a_p*Tp[inertial_deg] + q*Tp[2*inertial_deg]
            L_P.append(f)
    else:
        for P, a_p in ap_list:
            inertial_deg = 2 if P.is_inert() else 1
            if P in bad_primes:
                # bad reduction
                f = 1 - twist(P.sage_ideal()) * a_p*Tp[inertial_deg]
            else:
                # good reduction
                q = P.norm()
                # note that the following will automatically make the polynomial into a degree 0 term
                # in the case that the twist introduces bad reduction
                f = 1 - twist(P.sage_ideal()) * a_p*Tp[inertial_deg] + twist(P.sage_ideal())**2 * q*Tp[2*inertial_deg]
            L_P.append(f)

    # Use the local factors of the L-series to compute the Dirichlet
    # series coefficients of prime-power index.
    #coefficients = [0,1] + [0]*(bound-1)
    i = 0
    while i < len(ap_list):
        P = ap_list[i][0]
        if P.is_split():
            s = L_P[i] * L_P[i+1]
            i += 2
        else:
            s = L_P[i]
            i += 1
        p = P.p
        # We need enough terms t so that p^t > bound
        accuracy_p = int(math.floor(math.log(bound)/math.log(p))) + 1
        series_p = s.add_bigoh(accuracy_p)**(-1)
        for j in range(1, accuracy_p):
            an_list[p**j] = series_p[j]

    extend_multiplicatively(an_list)
    normalized_an_list = []
    for n, an in enumerate(an_list):
        if n == 0:
            continue
        else:
            normalized_an_list.append(an/RR(n).sqrt())

    #Q = RR(level.norm() * 25/(pi**2)).sqrt()
    Q = RR(level.norm()).sqrt() * 5/(4*pi**2)
    if twist is not None:
        Q = Q * twist.conductor().norm()
    poles = []
    residues = []
    period = 0


    L = Lfunction_D("", 2, normalized_an_list, period, Q, sign, [1,1], [.5,.5], poles, residues)
    return L.value(.5)

    

def make_tot_pos(n):
    """
    Given an element n of Q(sqrt5), return a totally positive unit multiple
    of n; i.e. return m such that m/n is a unit and m is taken to a positive
    number under both embeddings of K into R.

    EXAMPLES::
        
        sage: from psage.modform.hilbert.sqrt5.ellcurve import K, a, phi, phi2, make_tot_pos
        sage: N = 15 + 30 * a
        sage: M = make_tot_pos(N); M
        45*a + 30
        sage: M/N
        a
        sage: phi(M) > 0
        True
        sage: phi2(M) > 0
        True
        sage: K.ideal(N) == K.ideal(M)
        True
    """
    if phi(n) > 0:
        if phi2(n) < 0:
            return n * a
        else:
            return n
    else:
        if phi2(n) < 0:
            return n * a * (1 - a)
        else:
            return n * (1 - a)

def mixed_periods_from_curve(E):
    """
    Given an elliptic curve over Q(sqrt5) return the mixed periods of E and
    its Galois conjugate. This function exists mainly for testing purposes.

    EXAMPLES::

        sage: from psage.modform.hilbert.sqrt5.ellcurve import K, a, mixed_periods_from_curve
        sage: E = EllipticCurve(K, [0, -a, a, 0, 0])
        sage: mm, mp, pm, pp = mixed_periods_from_curve(E)
        sage: mm
        -8.10950089686...
        sage: mp
        15.7167643544...*I
        sage: pm
        23.8695845386...*I
        sage: pp
        46.260878468...
    """
    from sage.all import I, imag

    L1 = E.period_lattice(phi)
    L2 = E.period_lattice(phi2)

    if L1.real_flag == 1:
        m1 = imag(L1.basis()[1]) * I
    else:
        m1 = imag(L1.basis()[1] * 2)*I

    p1 = L1.basis()[0]

    if L2.real_flag == 1:
        m2 = imag(L2.basis()[1]) * I
    else:
        m2 = imag(L2.basis()[1] * 2)*I
    
    p2 = L2.basis()[0]

    return m1 * m2, m1 * p2, p1 * m2, p1 * p2
    
def find_rational(x, y, eps):
    """
    return the rational number a/b with smallest denominator
    such that real(x/y) = a/b up to an error of epsilon/b

    This is a stupid function. There must be a good way to do this.

    EXAMPLES::

        sage: from psage.modform.hilbert.sqrt5.ellcurve import find_rational
        sage: find_rational(.777, 1, .03)
        7/9
    """
    from sage.all import real

    X = real(x/y)
    b = Integer(1)
    while(  abs(X * b - round(X * b)) > eps ):
        b = b + 1
    return Integer(round(X*b))/b

def find_numerator(z):
    """
    Try to recognize z as a rational number, and return its numerator. Implementation: look at
    the convergents of the continued fraction expansion for z, and choose one just before
    the denominator jumps by a factor of more than 25, or the last one if we run out of precision

    EXAMPLES::
        
        sage: from psage.modform.hilbert.sqrt5.ellcurve import find_numerator
        sage: find_numerator(float(7/9))
        7
        sage: find_numerator(.22)
        2
    
    Note that our heuristic fails sometimes when there is a fairly good approximation
    to z with a smaller denomator than we are looking for.

        sage: find_numerator(float(4115/2263))
        20
        sage: float(4115/2263) - float(20/11)
        0.0002...

    What happens in this case is that there is a big jump in denominators in the continued
    fraction expansion:
        
        sage: continued_fraction(float(4115/2263))
        [1, 1, 4, 1, 1, 40, 1, 1, 1, 1]

    We see that 40, which causes us to stop.
    """

    C = continued_fraction(z)
    possibilities = C.convergents()
    # don't know the best way to do this. for now we are just going
    # to look through the list and stop if we hit the end or or the denominator
    # jumps by more than a factor of 25
    
    last_numerator = 0
    last_denominator = 1
    for x in possibilities:
        if x.denominator() > last_denominator * 25 and last_numerator != 0:
            return last_numerator
        else:
            last_denominator = x.denominator()
            last_numerator = x.numerator()

    return last_numerator

def recognize_rational(z):
    """
    Try to recognize z as a rational number, and return this rational. Implementation: look at
    the convergents of the continued fraction expansion for z, and choose one just before
    the denominator jumps by a factor of more than 25, or the last one if we run out of precision

    EXAMPLES::
        
        sage: from psage.modform.hilbert.sqrt5.ellcurve import find_numerator
        sage: find_numerator(float(7/9))
        7
        sage: find_numerator(.22)
        2
    
    Note that our heuristic fails sometimes when there is a fairly good approximation
    to z with a smaller denomator than we are looking for.

        sage: find_numerator(float(4115/2263))
        20
        sage: float(4115/2263) - float(20/11)
        0.0002...

    What happens in this case is that there is a big jump in denominators in the continued
    fraction expansion:
        
        sage: continued_fraction(float(4115/2263))
        [1, 1, 4, 1, 1, 40, 1, 1, 1, 1]

    We see that 40, which causes us to stop.
    """

    C = continued_fraction(z)
    possibilities = C.convergents()
    print possibilities
    # don't know the best way to do this. for now we are just going
    # to look through the list and stop if we hit the end or or the denominator
    # jumps by more than a factor of 25
    
    last_numerator = 0
    last_denominator = 1
    for x in possibilities:
        if x.denominator() > last_denominator * 25 and last_numerator != 0:
            return last_numerator/last_denominator
        else:
            last_denominator = x.denominator()
            last_numerator = x.numerator()

    return x

def compute_period_guesses(level, aplist, aplist_size, sign, ncpus = 0, limit1 = 5, chi_norm_bound = 1000, verbose=0):
    """
    Given some L-function data for an elliptic curve/hilbert modular form over Q(sqrt 5),
    multiples of its mixed periods.
    """

    I = CC(0,1)

    # we start by finding some characters that will suit out purposes.
    chi_mm_list, chi_mp_list, chi_pm_list, chi_pp_list = find_good_characters(sign, level, chi_norm_bound, limit=limit1)


    # and now with these characters we compute multiples of the mixed the periods of the L-functions

    primes = primes_of_bounded_norm(aplist_size)

    mm_list, mp_list, pm_list, pp_list = many_mixed_twisted_periods(zip(primes, aplist), [chi_mm_list, chi_mp_list, chi_pm_list, chi_pp_list], precision=53, ncpus=ncpus, verbose=verbose, level=level, sign=1)
    # note that the precision argument is currently unused, but we put it there anyway

    # the function for computing mixed periods is careless about a +- sign,
    # so for each nonzero multiple we adjust to get the correct sign and to strip the
    # zero real or imaginary part.
    def assign_sign(z, sign):
        # given a number z, adjust it's sign so that it
        # is the same as sign
        if z * sign > 0:
            return z
        else:
            return -z

    mm_list = [assign_sign(real(z), -1) for z in mm_list if z != 0]
    mp_list = [assign_sign(imag(z), 1) * I for z in mp_list if z != 0]
    pm_list = [assign_sign(imag(z), 1) * I for z in pm_list if z != 0]
    pp_list = [assign_sign(real(z), 1) for z in pp_list if z != 0]

    if len(mm_list) == 0 or len(mp_list) == 0 or len(pm_list) == 0 or len(pp_list) == 0:
        raise ValueError("""We couldn't find enough good characters with nonzero twists. If the level is not a square, try increasing limit1 and/or the chi_norm_bound.""")
  
    mm_ratio_list = [recognize_rational( (mm_list[0]/z).real() ) for z in mm_list]
    mp_ratio_list = [recognize_rational( (mp_list[0]/z).real() ) for z in mp_list]
    pm_ratio_list = [recognize_rational( (pm_list[0]/z).real() ) for z in pm_list]
    pp_ratio_list = [recognize_rational( (pp_list[0]/z).real() ) for z in pp_list]

    if verbose:
        print "Guesses for ratios of the multiples of mixed periods:"
        print "(If any of these numbers looks too big, something wrong probably happened.)"
        print "       ", mm_ratio_list
        print "       ", mp_ratio_list
        print "       ", pm_ratio_list
        print "       ", pp_ratio_list
        print "actual numbers:"
        print "       ", [mm_list[0]/z for z in mm_list]
        print "       ", [mp_list[0]/z for z in mp_list]
        print "       ", [pm_list[0]/z for z in pm_list]
        print "       ", [pp_list[0]/z for z in pp_list]
    
    mm_divider_guess = lcm([z.numerator() for z in mm_ratio_list])
    mp_divider_guess = lcm([z.numerator() for z in mp_ratio_list])
    pm_divider_guess = lcm([z.numerator() for z in pm_ratio_list])
    pp_divider_guess = lcm([z.numerator() for z in pp_ratio_list])

    mm = mm_list[0]/mm_divider_guess
    mp = mp_list[0]/mp_divider_guess
    pm = pm_list[0]/pm_divider_guess
    pp = pp_list[0]/pp_divider_guess

    if verbose:
        print "Computed the following period guesses:"
        print "   mm =", mm, "; mm1 =", mm_list[0]
        print "   mp =", mp, "; mp1 =", mp_list[0]
        print "   pm =", pm, "; pm1 =", pm_list[0]
        print "   pp =", pp, "; pp1 =", pp_list[0]

    return mm, mp, pm, pp

def find_curve_from_period_guesses(mm, mp, pm, pp, level, aplist, ncpus = 0, limit2 = (0,5), limit3 = (0,5), limit4 = 0, verbose=0, implementation = 'cython'):
    """
    Given multiples of the mixed periods of an elliptic curve over Q(sqrt 5), try to find that curve,
    and return it if found.
    """
    import sys
    I = CC(0, 1)

    level_gen = level.gens_reduced()[0]

    # if we knew the actual periods of the curve and the actual discriminant,
    # finding the curve would be no problem. We only know multiples of the periods,
    # however, and we only know the factors of the discriminant, so we will need
    # to iterate through the various possibilities to find out curve.

    D_base = level_gen

    D_base_embeddings = [phi(D_base), phi2(D_base)]
    D_factors = []
    for p, e in D_base.factor():
        if not 1/p in OK:
            D_factors.append(p)

    D_factors_embeddings = [ (phi(p), phi2(p)) for p in D_factors]
 
    # the following is all numbers a^n (1-a)^m for 0 <= n,m < 12. There are 44 of them,
    units = [-55*a + 89, 1, -55*a - 34, -21*a + 34, 3*a - 5, a - 1, 3*a + 2, -13*a - 8, 5*a - 8, -34*a + 55, 21*a + 13, -a - 1, 5*a + 3, a - 2, a, a + 1, 34*a + 21, -13*a + 21, 2*a - 3, 2*a + 1, 13*a - 21, -a, -a + 1, -a + 2, -2*a + 3, 89*a + 55, -34*a - 21, 8*a - 13, -5*a - 3, -2*a - 1, -21*a - 13, -5*a + 8, 13*a + 8, 34*a - 55, 8*a + 5, -3*a - 2, -8*a - 5, -3*a + 5, 21*a - 34, -89*a + 144, 55*a + 34, 55*a - 89, -1, -8*a + 13]

    units_embeddings = [(phi(u), phi2(u)) for u in units]

    # for now we are going to ignore the -- period, since we don't actually use it for anything.
    # it could be better to use it to try to find a different guess for the ++ period, though.

    D_guesses = []
    for N in range(*limit2):
        for exponent_guesses in IntegerVectors(N, len(D_factors)):
            for u in units:
                D_guesses.append( u * D_base * prod( [p**e for (p,e) in zip(D_factors, exponent_guesses)] ) )

    if verbose > 1:
        print "Discriminant guesses:"
        for D in D_guesses:
            print "   ", D

    D_guess_embeddings = [ (phi(x), phi2(x)) for x in D_guesses ]

    if verbose:
        print "built list of", len(D_guesses), "D guesses."
        sys.stdout.flush()

    # j1 will be phi(j(E)) and j2 will be phi2(j(E))

    j1a_possibilities = []
    j1b_possibilities = []
    j2a_possibilities = []
    j2b_possibilities = []

    J = j_invariant_qexp(30)

    C = ComplexField(200)
    for n in range(*limit3):
        for j, k, l in IntegerVectors(n, 3):
            j += 1
            k += 1
            l += 1
            if gcd( [j, k, l] ) > 1:
                continue

            tau1a = reduce_tau( mp/pp * l/j )[0]
            tau1b = reduce_tau( .5*(1 + mp/pp * l/j ))[0]
            tau2a = reduce_tau( pm/pp * l/k )[0]
            tau2b = reduce_tau( .5*(1 + pm/pp * l/k ))[0]

            j1a_possibilities.append( J( C(2 * pi * I * tau1a).exp().real() ) )
            j1b_possibilities.append( J( C(2 * pi * I * tau1b).exp().real() ) )
            j2a_possibilities.append( J( C(2 * pi * I * tau2a).exp().real() ) )
            j2b_possibilities.append( J( C(2 * pi * I * tau2b).exp().real() ) )

    if verbose:
        print "built list of", len(j1a_possibilities), "j invariant guesses."
        sys.stdout.flush()

    if verbose > 1:
        print "J possibilities:"
        for x, y, z, w in zip(j1a_possibilities, j1b_possibilities, j2a_possibilities, j2b_possibilities):
            print "   ", x, y, z, w

    # with our list of discriminant guesses j-invariant guesses we now proceed to try to find
    # the curve

    sqrt5 = RR(sqrt(5))

    total_guesses = len(j1a_possibilities) * len(D_guesses)
    count = 0
    
    if implementation=="cython":
        E = search_possibilities(D_guesses, D_guess_embeddings, j1a_possibilities, j1b_possibilities, j2a_possibilities, j2b_possibilities, level, aplist, verbose)
        return E

    for n in range(len(D_guesses)):
        D = D_guesses[n]
        D1, D2 = D_guess_embeddings[n]
        if D1 > 0:
            j1_possibilities = j1a_possibilities
        else:
            j1_possibilities = j1b_possibilities
        if D2 > 0:
            j2_possibilities = j2a_possibilities
        else:
            j2_possibilities = j2b_possibilities 
        for j1, j2 in zip(j1_possibilities, j2_possibilities):
            if verbose:
                if count % 100 == 0:
                    print "Completed", count, "guesses out of", total_guesses
            count += 1
            c4a = (j1 * D1).nth_root(3)
            c4b = (j2 * D2).nth_root(3)

            if verbose > 2:
                print "Trying c4 = "
                print "  phi(c4) =", c4a
                print " phi2(c4) =", c4b
                sys.stdout.flush()

            # now write c4 as A + B sqrt(5)
            A = ZZ(round(c4a + c4b))/2
            B = ZZ(round((c4a - c4b)/sqrt5))/2

            #print "trying c4 around", A, "+", B,"sqrt5"
            base_x = ZZ(round(A - B)) # A - B should be an integer, but
                              # our precision isn't perfect, so it might not be
            base_y = ZZ(2 * B)
            if verbose > 2:
                print "trying c4 around", A, "+", B,"sqrt5 =", base_x + base_y * a
            #c4_possibilities = [x + a * y for (x,y) in cartesian_product_iterator([srange(base_x - 5, base_x + 5), srange(base_y - 5, base_y + 5)])]

            if limit4 == 0:
                c4_possibilities = [base_x + base_y * a]
            else:
                c4_possibilities = [x + a * y for (x,y) in cartesian_product_iterator([srange(base_x - limit4, base_x + limit4), srange(base_y - limit4, base_y + limit4)])]
                

            for c4 in c4_possibilities:
                #try:
                #    _ = ZZ(c4)
                #except TypeError:
                #print "Trying c4 =", c4
                Z = c4**3 - 1728 * D
                if Z.is_square():
                    c6 = sqrt(Z)
                    if verbose > 1:
                        print "Finding c6=", c6
                    E = EllipticCurve_from_c4c6(c4, c6)
                    if verbose > 1:
                        print
                    E = E.integral_model()
                    if verbose > 1:
                        print "Found a curve!!!!", c4, c6, E.conductor()
                    if E.conductor() == level:
                        E = E.global_minimal_model()
                        z = compute_aplist(E, 100)
                        if z == aplist[:len(z)]:
                            return E

                    c6 = -c6
                    E = EllipticCurve_from_c4c6(c4, c6)
                    if verbose > 1:
                        print
                    E = E.integral_model()
                    if verbose > 1:
                        print "Found a curve!!!!", c4, c6, E.conductor()
                    if E.conductor() == level:
                        E = E.global_minimal_model()
                        z = compute_aplist(E, 100)
                        if z == aplist[:len(z)]:
                            return E
                        #return E
                    if verbose > 1:
                        print
                else:
                    #print "that c4 is no good"
                    pass

         



def find_curve_from_Lfunction(level, aplist, aplist_size, sign, ncpus = 1, limit1=1, limit2 = (0,5), limit3 = (0,5), limit4 = 0, verbose=0, chi_norm_bound = 1000, implementation = "cython"):
    """
    Given some $L$-function data, find an elliptic curve over $\Q(\sqrt 5)$
    with this $L$-function. Uses the method described in Lassina Dembele's
    paper "..."

    This method has not really been tested enough yet to know what the best
    way to use it is. In general, if we know a lot ap values, then it should
    be a good thing to set limit1 to be something like 4 or 5, or maybe even
    higher. However, if we don't know that many coefficients, setting limit1
    this high will causes problems because of precision issues. The reason
    for setting limit1 higher, though, is so that limit3 can be set lower, so
    that a much smaller space needs to be searched.

    In principal, one should be able to just so something like

    E = find_curve_from_Lfunction(level, aplist, aplist_size, sign)

    (with those variables set appropriately, of course.) If it doesn't work,
    you can try again with limit2 and/or limit3 set to a larger range. If
    the ranges are large enough and mixed periods were computed to enough
    precision, then the curve should be found. If there isn't enough
    precision, then it might be necessary to try increasing limit4. (Or
    you might roll some dice and use those to pick coefficients
    of an elliptic curve. Maybe it will be the one you are looking for!)

    In practice, it is probably best to take limit1 to be something between
    3 and 10, set verbose = 1, and look at the output to make sure that
    nothing bad happened. If the code is having trouble automatically
    identifying rational numbers because of limited precision, you can look
    at the output and try to do it manually, and then use the function
    find_curve_from_period_guesses(). Or you can just set limit1 to be lower
    and set limit3 higher, and let the computer do the work.

    INPUT::
        level - An ideal. The conductor of the curve/the level of the Hilbert
            modular form. Must not be a square.
        
        aplist - A list of Fourier coefficients of the modular form at primes.
            The nth entry in this list should correspond to the nth entry in
            primes_of_bounded_norm.

        aplist_size - The norm of prime corresponding the the last entry in aplist

        sign - The sign of the functional equation

        chi_norm_bound - the largest norm conductor we are willing to twist by

        ncpu - The number of cpus to use while computing twists of the L-function

        limit1 - the maximum number of twists we compute for each mixed period

        limit2 - a range on the powers of the discriminant we are going to
            try. We will try all possibilities D = level * prod(p_k^e_k) * u
            where p_k are the divisors of the level and sum e_k is in the range
            limit2

        limit3 - a range of the integers we will try dividing the mixed period
            guesses by. After computing initial guesses, we will try guesses
            of the form (\Omega_k)/m_k where \sum m_k is in the range limit3

        limit4 - only used in the sage implementation of searching after
            computing periods. After computing an initial guess of
            c4 = A + B \varphi, we will try all possibilities in the range
            (A += limit4), (B += limit4). This greatly slows things down
            (by a factor of 4 limit4^2, and probably shouldn't be necessary
            if we have enough precision

        implementation - 'sage' or 'cython'. The code we will use for
            searching the space of possible guesses. The cython implemenation
            is much faster, but has some overflow issues which may make it
            miss things that it should find. (But I don't know yet of an
            example where this happens. It is possible that when we have
            such large numbers we shouldn't be trying to find the curve
            anyway.)

    EXAMPLE::
        
        A simple example from a curve that we know about comes first.

        sage: from psage.modform.hilbert.sqrt5.ellcurve import *
        sage: E = EllipticCurve(K, [0, -a, a, 0, 0])
        sage: E2 = find_curve_from_Lfunction(E.conductor(), compute_aplist(E, 40000), 40000, 1, ncpus=4, limit1=3, verbose=0) # not really random output, but too much right now
        
        The curve we get back is in fact isomorphic to the original, but in general that is asking too much,
        so we just check isogeny.

        sage: E.is_isogenous(E2)
        True

        Here is a nice example that works well, which relies on a bunch of Fourier coefficients that
        William Stein computed for a Hilbert modular form of level 48 - 8*a.

        sage: aplist = load("http://wstein.org/home/wstein/reu/2011/aplist-1856b.sobj") # some random message about loading gets printed here
        sage: level = K.ideal(48 - 8*a)
        sage: sign = 1
        sage: E = find_curve_from_Lfunction(level, aplist, 29000, sign, ncpus=4, limit1=3, limit2=(0, 10), limit3=(0, 5), limit4 = 0, verbose=1) # again, not really random, but noisy
        sage: E
        Elliptic Curve defined by y^2 = x^3 + (a+1)*x^2 + (-1439*a-893)*x + (-32764*a-20247) over Number Field in a with defining polynomial x^2 - x - 1
        sage: E.conductor() == level
        True
        sage: x = compute_aplist(E, 100); x == aplist[:len(x)]
        True
    """
    
    
    # we start by guessing multiples of the mixed periods of the curve.
    mm, mp, pm, pp = compute_period_guesses(level, aplist, aplist_size, sign, ncpus, limit1, verbose=verbose, chi_norm_bound=chi_norm_bound)

    # and then we find the curve from those guesses
    E = find_curve_from_period_guesses(mm, mp, pm, pp, level, aplist, ncpus, limit2 = limit2, limit3 = limit3, limit4=limit4, verbose=verbose, implementation = implementation)

    if E is None:
        print "did not find curve, returning period guesses."
        return mm, mp, pm, pp
    # it is as simple as that!
    return E



def find_curve_from_curve2(E, precision=20, verbose = 0, print_period_multiples = False, ncpus = 0, limit1=3, limit2 = 10, number_of_coefficients = 40000):
    """
    Given an elliptic curve E over Q(sqrt5), we go through the steps of Dembele's algorithm,
    except that we use the information from the curve to avoid guessing.

    This means that we:
        1. Compute the global root number for the curve.
        2. Find good characters to twist by.
        3. Twist by a few of those characters to compute multiples of the mixed periods
           of the curve.
        4. Use ratios of differnet multiples of the same mixed period to discover
           which multiple of the mixed period the first one (probably) is.
        5. Use these periods and (cheat) the determinant of the curve to compute the
           j-invariant of the curve and its conjugate, and from the j-invariants the c4
           of the curve and its conjugate. This gives us c4, and we use the determinant
           to get c6.

        It should be possible to get around the guessing in steps 4 and 5 by trying many different
        values until we actually find a curve. These are the only places where we use any information
        in the curve that is not directly contained in its L-function. 

    EXAMPLE::
        sage: from psage.modform.hilbert.sqrt5.ellcurve import *
        sage: E = EllipticCurve(K, [0,0,a,27*a-43,-80*a+128])

    """
    from sage.all import I
    import sys
    mmE, mpE, pmE, ppE = mixed_periods_from_curve(E)

    level = E.conductor()
    disc = E.discriminant()

    N = level.norm()
    level_gen = make_tot_pos(level.gens_reduced()[0])

    start_time = time.time()

    #an_list = build_an_list(E, 1500)
    #f = {'an_list':an_list, 'level_gen':level_gen}
    if verbose:
        print "finding characters"
    sys.stdout.flush()

    L = lseries_dokchitser(E, prec=20) # we are just going to use this to determine the sign
                                       # of the functional equation (global root number)

    chi_mm_list, chi_mp_list, chi_pm_list, chi_pp_list = find_good_characters(L.eps, level, bound=1500, limit=limit1)

    if verbose:
        print "computing mixed periods"
    sys.stdout.flush()

    ap_list = compute_aplist(E, number_of_coefficients)
    primes = primes_of_bounded_norm(number_of_coefficients)
    ap_list = zip(primes, ap_list)

    mm_list, mp_list, pm_list, pp_list = many_mixed_twisted_periods(ap_list, [chi_mm_list, chi_mp_list, chi_pm_list, chi_pp_list], precision=precision, ncpus=ncpus, verbose=verbose, level=level, sign=1)

    if verbose:
        print "computed periods. time so far:", time.time() - start_time
        sys.stdout.flush()

    #mm1, mp1, pm1, pp1 = [_2 for (_1, _2) in sorted(list(mixed_period( [ (_, f, chi, 30000, 0) for (_, chi) in enumerate([chi_mm_list[0], chi_mp_list[0], chi_pm_list[0], chi_pp_list[0]])] )))]

    mm1 = 0
    mp1 = 0
    pm1 = 0
    pp1 = 0

    def grab_first_nonzero(M):
        z = 0
        while z == 0:
            print z
            z = M.pop(0)

        return z

    def assign_sign(z, sign):
        # given a number z, adjust it's sign so that it
        # is the same as sign
        if z * sign > 0:
            return z
        else:
            return -z

    mm_list = [assign_sign(real(z), -1) for z in mm_list if z != 0]
    mp_list = [assign_sign(imag(z), 1) * I for z in mp_list if z != 0]
    pm_list = [assign_sign(imag(z), 1) * I for z in pm_list if z != 0]
    pp_list = [assign_sign(real(z), 1) for z in pp_list if z != 0]

    mm1, mp1, pm1, pp1 = (grab_first_nonzero(z) for z in (mm_list, mp_list, pm_list, pp_list))

    print mm1, mp1, pm1, pp1
    print (mm_list, mp_list, pm_list, pp_list)
    print [find_numerator(mm1/z) for z in mm_list]
    print [find_numerator(mp1/z) for z in mp_list]
    print [find_numerator(pm1/z) for z in pm_list]
    print [find_numerator(pp1/z) for z in pp_list]
    
    mm_divider_guess = lcm([find_numerator(mm1/z) for z in mm_list])
    mp_divider_guess = lcm([find_numerator(mp1/z) for z in mp_list])
    pm_divider_guess = lcm([find_numerator(pm1/z) for z in pm_list])
    pp_divider_guess = lcm([find_numerator(pp1/z) for z in pp_list])

    mm = mm1/mm_divider_guess
    mp = mp1/mp_divider_guess
    pm = pm1/pm_divider_guess
    pp = pp1/pp_divider_guess


    if verbose:
        print "starting with periods:"
        print "   mm =", mm1, "/", mm_divider_guess, "=", mm, " ; actual:", mmE
        print "   mp =", mp1, "/", mp_divider_guess, "=", mp, " ; actual:", mpE
        print "   pm =", pm1, "/", pm_divider_guess, "=", pm, " ; actual:", pmE
        print "   pp =", pp1, "/", pp_divider_guess, "=", pp, " ; actual:", ppE
        print ""
        print "Actual j-invariants:"
        print "   phi(j(E)) =", phi(E.j_invariant())
        print "  phi2(j(E)) =", phi2(E.j_invariant())
        print "Actual c4:"
        print "        c4 =", E.c4()
        print "   phi(c4) =", phi(E.c4())
        print "  phi2(c4) =", phi2(E.c4())
        print "   minpoly =", E.c4().minpoly()
        sys.stdout.flush()


    level_factors = [p.gens_reduced()[0] for (p,e) in factor(level)]
    number_of_factors = len(level_factors)

    
    unit_multipliers = []

    for n in range(-6, 6):
        for m in range(-6, 6):
            unit_multipliers.append(a**n * (1-a)**m)

    unit_multipliers = list(Set(unit_multipliers))

    print unit_multipliers

    def disc_guess_generator():
        for n in range(50):
            for power_list in IntegerVectors(n, number_of_factors):
                D_guess = level_gen * prod((p**e for (p,e) in zip(level_factors, power_list)))
                for u in unit_multipliers:
                    yield D_guess * u
    
    # the following strategy for looping through guesses should probably be thought
    # about more carefully. we will do something like: take an integer vector with
    # ... (THINK SOME MORE ON THIS.)

    @parallel(ncpus=ncpus)
    def try_disc(D):
        #for u in unit_multipliers:
        u = 1
        for m in range(limit2):
            for _divider_list in IntegerVectors(m, 4):
                divider_list = [x + 1 for x in  _divider_list]
                if gcd(divider_list) == 1:
                    mm_try = mm/divider_list[0]
                    mp_try = mp/divider_list[1]
                    pm_try = pm/divider_list[2]
                    pp_try = pp/divider_list[3]
                    if abs(abs(mm_try * pp_try/(pm_try * mp_try)) - 1) < .001:
                        E_list = try_to_find_curve(mm_try, mp_try, pm_try, pp_try, u * D, level, verbose=verbose)
                        for E2 in E_list:
                            E2 = E2.integral_model().global_minimal_model()
                            if anlist(E, 100) == anlist(E2, 100):
                                return E2

    for _ in try_disc( disc_guess_generator() ):
        E2 = _[1]
        if E2 is not None:
            print "returning curve. elapsed time:", time.time() - start_time
            return E2

    #for n in range(50): #this is a very big range. don't expect to go that far
    #    for power_list in IntegerVectors(n, number_of_factors):
    #        D_guess = level_gen * prod((p**e for (p,e) in zip(level_factors, power_list)))
    #        for u in unit_multipliers:
    #            for m in range(limit2):
    #                for _divider_list in IntegerVectors(m, 4):
    #                    divider_list = [x + 1 for x in  _divider_list]
    #                    if gcd(divider_list) == 1:
    #                        mm_try = mm/divider_list[0]
    #                        mp_try = mp/divider_list[1]
    #                        pm_try = pm/divider_list[2]
    #                        pp_try = pp/divider_list[3]
    #                        if abs(abs(mm_try * pp_try/(pm_try * mp_try)) - 1) < .0001:
    #                            E2 = try_to_find_curve(mm_try, mp_try, pm_try, pp_try, D_guess * u, level, verbose=verbose)
    #                            if E2 is not None:
    #                                if verbose:
    #                                    print "returning curve. elapsed time:", time.time() - start_time
    #                                return E2
                                    
    if verbose:
        print "no curve found. elapsed time:", time.time() - start_time


def find_curve_from_curve(E, precision=20, verbose = 0, print_period_multiples = False):
    """
    This function is rendered obsolete by later developments.  Eventually the same
    function with a number 2 will replace this.


    Given an elliptic curve E over Q(sqrt5), we go through the steps of Dembele's algorithm,
    except that we use the information from the curve to avoid guessing.

    This means that we:
        1. Compute the global root number for the curve.
        2. Find good characters to twist by.
        3. Twist by those characters to compute multiples of the mixed periods
           of the curve.
        4. (Cheat) Compute the actual mixed periods of the curve to determine
           which multiples of the periods we've computed.
        5. Use these periods and (cheat) the determinant of the curve to compute the
           j-invariant of the curve and its conjugate, and from the j-invariants the c4
           of the curve and its conjugate. This gives us c4, and we use the determinant
           to get c6.

        It should be possible to get around the guessing in steps 4 and 5 by trying many different
        values until we actually find a curve. These are the only places where we use any information
        in the curve that is not directly contained in its L-function. 

    """
    from sage.all import I
    import sys
    mmE, mpE, pmE, ppE = mixed_periods_from_curve(E)

    level = E.conductor()
    disc = E.discriminant()

    N = level.norm()
    level_gen = make_tot_pos(level.gens_reduced()[0])


    #an_list = build_an_list(E, 1500)
    #f = {'an_list':an_list, 'level_gen':level_gen}
    if verbose:
        print "finding characters"
    sys.stdout.flush()

    L = lseries_dokchitser(E, prec=15) # we are just going to use this to determine the sign
                                       # of the functional equation (global root number)

    chi_mm_list, chi_mp_list, chi_pm_list, chi_pp_list = find_good_characters(L.eps, level, bound=1500, limit=10)

    if verbose:
        print "computing mixed periods"
    sys.stdout.flush()



    #mm1, mp1, pm1, pp1 = [_2 for (_1, _2) in sorted(list(mixed_period( [ (_, f, chi, 30000, 0) for (_, chi) in enumerate([chi_mm_list[0], chi_mp_list[0], chi_pm_list[0], chi_pp_list[0]])] )))]

    mm1 = 0
    mp1 = 0
    pm1 = 0
    pp1 = 0

    mm_n = 0
    mp_n = 0
    pm_n = 0
    pp_n = 0
    while(abs(mm1) < .1):
        if verbose:
            print "computing mixed twisted period with chi.conductor().norm =", chi_mm_list[mm_n].conductor().norm()
            sys.stdout.flush()
        mm1 = mixed_twisted_period(E, chi_mm_list[mm_n], precision=precision)
        if verbose:
            print mm1
        #mm1 = mixed_periods_from_curve (1, f, chi_mm_list[mm_n], 30000, 0)
        mm_n = mm_n + 1

    while(abs(mp1) < .1):
        if verbose:
            print "computing mixed twisted period with chi.conductor().norm =", chi_mp_list[mp_n].conductor().norm()
            sys.stdout.flush()
        mp1 = mixed_twisted_period(E, chi_mp_list[mp_n], precision=precision)
        if verbose:
            print mp1
        #mp1 = mixed_period(1, f, chi_mp_list[mp_n], 30000, 0)
        mp_n = mp_n + 1

    while(abs(pm1) < .1):
        if verbose:
            print "computing mixed twisted period with chi.conductor().norm =", chi_pm_list[pm_n].conductor().norm()
            sys.stdout.flush()
        pm1 = mixed_twisted_period(E, chi_pm_list[pm_n], precision=precision)
        if verbose:
            print pm1
        #pm1 = mixed_period(1, f, chi_pm_list[pm_n], 30000, 0)
        pm_n = pm_n + 1

    while(abs(pp1) < .1):
        if verbose:
            print "computing mixed twisted period with chi.conductor().norm =", chi_pp_list[pp_n].conductor().norm()
            sys.stdout.flush()
        pp1 = mixed_twisted_period(E, chi_pp_list[pp_n], precision=precision)
        if verbose:
            print pp1
        #pp1 = mixed_period(1, f, chi_pp_list[pp_n], 30000, 0)
        pp_n = pp_n + 1

    mm1 = real(mm1)
    mp1 = imag(mp1) * I
    pm1 = imag(pm1) * I
    pp1 = real(pp1)

    if mm1 > 0:
        mm1 = -mm1
    if imag(mp1) < 0:
        mp1 = -mp1
    if imag(pm1) < 0:
        pm1 = -pm1
    if pp1 < 0:
        pp1 = -pp1

    if verbose:
        print "computed mixed periods: "
        print "  mm1 =", mm1, "; conductor norm", chi_mm_list[0].conductor().norm()
        print "  mp1 =", mp1, "; conductor norm", chi_mp_list[0].conductor().norm()
        print "  pm1 =", pm1, "; conductor norm", chi_pm_list[0].conductor().norm()
        print "  pp1 =", pp1, "; conductor norm", chi_pp_list[0].conductor().norm()
        print ""
        print "actual periods:"
        print "   mm =", mmE
        print "   mp =", mpE
        print "   pm =", pmE
        print "   pp =", ppE
        print "finding which multiples we've found"
    sys.stdout.flush()

    mm_divider = find_rational(mm1, mmE, .05)
    mp_divider = find_rational(mp1, mpE, .05)
    pm_divider = find_rational(pm1, pmE, .05)
    pp_divider = find_rational(pp1, ppE, .05)

    mm = mm1/mm_divider
    mp = mp1/mp_divider
    pm = pm1/pm_divider
    pp = pp1/pp_divider

    if print_period_multiples:
        print mm_divider, mp_divider, pm_divider, pp_divider

    if verbose:
        print "using periods:"
        print "   mm =", mm1, "/", mm_divider, "=", mm
        print "   mp =", mp1, "/", mp_divider, "=", mp
        print "   pm =", pm1, "/", pm_divider, "=", pm
        print "   pp =", pp1, "/", pp_divider, "=", pp
        print ""
        print "Actual j-invariants:"
        print "   phi(j(E)) =", phi(E.j_invariant())
        print "  phi2(j(E)) =", phi2(E.j_invariant())
        print "Actual c4:"
        print "        c4 =", E.c4()
        print "   phi(c4) =", phi(E.c4())
        print "  phi2(c4) =", phi2(E.c4())
        print "   minpoly =", E.c4().minpoly()

    E2 = try_to_find_curve(mm, mp, pm, pp, [disc], level, verbose=verbose)
    return E2


def build_an_list(E, bound = 2000, verbose=0):
    import sys
    level = E.conductor()

    an_list = {}
    ap_dict = {}

    #print "finding ap values and creating prime ideals"
    sys.stdout.flush()

    ideals = K.ideals_of_bdd_norm(bound)
    prime_ideals = primes_of_bounded_norm(bound)
    ap_list = compute_aplist(E, bound)
    for prime, ap in zip(prime_ideals, ap_list):
        ap_dict[prime.sage_ideal()] = ap
    #for n in range(1, bound):
    #    if verbose:
    #        print "looking for primes of norm", n
    #    for p in ideals[n]:
    #        p = K.ideal(p.gens_reduced())
    #        #if is_prime(p):
    #        #    ap_dict[p] = psage.ellcurve.lseries.aplist.ap(E, p)
    #        if is_prime(p):
    #            if p.divides(E.discriminant()):
    #                #print "computing ap for p =", p, "n =", n
    #                ap_dict[p] = psage.ellcurve.lseries.aplist.ap(E, p)
    #            else:
    #                #print "computing ap for p =", p, "n =", n
    #                ap_dict[p] = E.change_ring(p.residue_field()).trace_of_frobenius()
                
    #print "finding an values"
    sys.stdout.flush()

    if 0: # bug somewhere in this code...
        for n in range(1, bound):
            for I in ideals[n]:
                an_list[I] = 1


        for p, ap in zip(prime_ideals, ap_list):
            p = p.sage_ideal()
            pk = 1
            for k in range(1, floor(log(bound)/log(p.norm()))):
                pk = pk * p
                if p.divides(level):
                    apk = ap**k
                else:
                    if k == 1:
                        apk = ap
                    else:
                        apk = ap * an_list[pk/p] - p.norm() * an_list[pk/(p * p)]
                for n in range(1, floor(bound/pk.norm())):
                    for I in ideals[n]:
                        if not p.divides(I):
                            print (pk*I).norm()
                            an_list[pk * I] *= apk

        # finally, since we want our an list to be indexed by totally positive generators,
        # and in slightly different form, we put it in that form

        an_dict = {}
        for n, an in an_list.iteritems():
            N = make_tot_pos(n.gens_reduced()[0])
            an_dict[N] = (phi(N), phi2(N), an)

        return an_dict
    else:
        for n in range(1, bound):
            if verbose:
                print "computing an for norm", n
            for I in ideals[n]:
                #print "computing an for norm", n
                an = 1
                F = I.factor()
                for (p, e) in F:
                    p = K.ideal(p.gens_reduced())
                    ap = ap_dict[p]
                    if p.divides(level):
                        an *= ap**e
                    else:
                        if e == 1:
                            an *= ap
                        elif e == 2:
                            an *= (ap**2 - p.norm())
                        elif e == 3:
                            an *= (ap**3 - 2 * p.norm() * ap)

                        elif e == 4:
                            an *= (ap**4 - 4 * ap**2 * p.norm() + p.norm()**2)
                        elif e == 5:
                            an *= ap**5 - 4 * ap**3 * p.norm() + 3 * ap * p.norm()**2
                        elif e == 6:
                            an *= ap**6 - 5*ap**4*p.norm() + 6*ap**2*p.norm()**2 - p.norm()**3
                        elif e == 7:
                            an *= ap**7 - 6*ap**5*p.norm() + 10*ap**3*p.norm()**2 - 4*ap*p.norm()**3
                        else:
                            raise ValueError("not ready to handle thing that are divisible by higher than cubic powers")
                N = make_tot_pos(I.gens_reduced()[0])
                an_list[N] = (phi(N), phi2(N), an)
        return an_list



def find_good_characters(sign, level, bound=100, limit=5):
    mm_list = []
    mp_list = []
    pm_list = []
    pp_list = []
    
    level_gen = make_tot_pos(level.gens_reduced()[0])

    for p in prime_range(bound):
        #print "looking for a character mod", p
        if p == 5 or p == 2:
            continue
        P_list = K.primes_above(p)
        N = P_list[0].norm()
        if N > bound + 100:
            continue
        for P in P_list:
            chi = QuadraticCharacter(P)
            if chi(-level_gen) != sign:
                continue
            if chi(a) == 1 and chi(1-a) == 1:
                pp_list.append((N,chi))
            elif chi(a) == 1 and chi(1-a) == -1:
                mp_list.append((N, chi))
            elif chi(a) == -1 and chi(1-a) == 1:
                pm_list.append((N, chi))
            else:
                mm_list.append((N, chi))

    pp_list = [w for (_, w) in sorted(pp_list)]
    pm_list = [w for (_, w) in sorted(pm_list)]
    mp_list = [w for (_, w) in sorted(mp_list)]
    mm_list = [w for (_, w) in sorted(mm_list)]
    
    return mm_list[:limit], mp_list[:limit], pm_list[:limit], pp_list[:limit]

def mixed_twisted_period(E, chi, precision = 20, verbose=0, level = None, sign = None):
    """
    return L(E, chi, 1) multiplied by some factors so that it is an integer (or rational?) multiple
    of a mixed period of E

    TODO: come up with a better name for this function
    """
    
    if isinstance(E, (list, tuple)):
        # assume that E is a list of ap values
        Lone = Lone_twist(E, sign, level, chi)
    else:
        try:
            L = lseries_dokchitser(E, twist=chi, prec=precision)
            Lone = L(1)
            if verbose:
                print "Computed L(1) to %d bits of precision using %d coefficients." % (precision, L.num_coeffs())
        except RuntimeError:
            # we need to be smarter here. the precision need to be either increased or
            # decreased according to what raised the RuntimeError
            Lone = 0
    

    # We only get the gauss sum correct up to +- 1, but that is an integer
    # Also, we will be able to correct the sign easily later if it is wrong.
    gauss_sum_term = RR(chi.conductor().norm()).sqrt()
    if chi(-1) == -1:
        gauss_sum_term *= CC(0, 1)


    return_value =  Lone * CC(sqrt(5) * gauss_sum_term) # TODO: figure out why we aren't multiplying by
                                               # 4pi^2. I mean, I never understood why we should
                                               # multiply by 4pi^2, but it worked with Lassina's
                                               # formula for L(1, f). So either his definition/formula
                                               # is wrong, or there is a difference of 4pi^2
                                               # between the definition of the modular form
                                               # L-function and the Elliptic Curve L-function
    
    tiny = .01 # return 0 if it seems like the L function is 0. This could be tricky,
               # because lcalc gives a rather small amount of precision with high conductor
               # and only 40000 coefficients, which is what our target is. In pratice the nonzero
               # numbers don't seem to be anywhere close to 0, though.
    if abs(return_value) < tiny:
        return 0
    else:
        return return_value

    #return L(1) * CC(-4 * pi**2 * sqrt(5) * gauss_sum_term)

def many_mixed_twisted_periods(E, chi_list, precision = 20, ncpus = 0, verbose=0, level=None, sign=None):
    """
    return L(E, chi, 1) multiplied by some factors so that it is an integer (or rational?) multiple
    of a mixed period of E

    For convenience, if chi_list is a list of lists, will return a list of lists, with each
    entry properly corrensponding to the input

    TODO: come up with a better name for this function
    """

    @parallel(ncpus=ncpus)
    def f(n):
        return mixed_twisted_period(E, chi_list[n], precision=precision, verbose=verbose, level=level, sign=sign)
    
    if isinstance(chi_list[0], (list, tuple)):
        list_lengths = [len(x) for x in chi_list]
        indices = [0]
        for n in list_lengths:
            indices.append(indices[-1] + n)
        concatenated_list = []
        for x in chi_list:
            concatenated_list.extend(x)

        chi_list = concatenated_list
        periods = [period for (_, period) in sorted(list(f(range(len(chi_list)))))]
        answer = []
        for previous, current in zip(indices, indices[1:]):
            answer.append(periods[previous:current])

    else:
        answer = [period for (_, period) in sorted(list(f(range(len(chi_list)))))]
    return answer


@parallel(ncpus=8)
def mixed_period(blah, f, chi, normbound = 10000000, verbose = 0):
    an_list = f['an_list']
    level_gen = f['level_gen']
    # right now f is the form with the above ap list
    if verbose:
        print "mixed_period() called with chi.conductor.norm() ==", chi.conductor.norm()

    N = phi(level_gen)      # a totally positive generator for the the level of f
    Nbar = phi2(level_gen)
    

    c = make_tot_pos(chi.conductor().gens_reduced()[0])
    nu = phi(c)
    nubar = phi2(c)
    
    D = RR(5)
    
    S = 0
    
    x,y,z,w = SR.var('x,y,z,w')    
    A = fast_callable( 1 - (-2 * pi * x * y/(nu * sqrt(D * N)).exp()), vars = [x,y], domain = RR)
    B = fast_callable( ( 2 * pi / sqrt(D) * (x * y / (nubar * sqrt(Nbar)) - z * w/(nu * sqrt(N))) ).exp(), vars = [x,y,z,w], domain = RR)
    
    for (n, (n1, n2, an)) in an_list.iteritems():
        if n.norm() > normbound:
            continue
        #n = make_tot_pos(n.gens_reduced()[0])
        S2 = 0

        #A = fast_callable( 1 - exp(-2 * pi * phi(n) * x/(nu * sqrt(D * N))), vars = [x], domain = R)
        #B = fast_callable( exp( 2 * pi / sqrt(D) * (x * y / (nubar * sqrt(Nbar)) - z * w/(nu * sqrt(N))) ), vars = [x,y,z,w], domain = R)
        for k in range(-20, 20): # not sure what kind of range to use on k right now
            #A = R(-2 * pi * phi(n) * phi(a)^(2 * k + 1) / (nu * sqrt(D * N)))
            #B = R(2 * pi / sqrt(D) * (phi2(n) * phi(1 - a)^(2 * k + 1) / (nubar * sqrt(Nbar)) - phi(n) * phi(a)^(2 * k)/(nu * sqrt(N))))
            S2 = S2 + A(n1, a1**(2 * k + 1))*B(n2, a2**(2 * k + 1), n1, a1**(2 * k))

        S = S + an * chi(n)/n.norm() * S2

    S = -S * 5/RR(2 * pi**2)
    
    return CC(-4 * pi**2 * sqrt(5) * chi.gauss_sum()) * S/5


def try_to_find_curve(mm, mp, pm, pp, D, level, verbose = 0):
    from sage.all import i

    # make sure everything has the correct sign:
    if mm > 0:
        mm = -mm
    if imag(mp) < 0:
        mp = -mp
    if imag(pm) < 0:
        pm = -pm
    if pp < 0:
        pp = -pp

    if verbose:
        print "Trying D =", D
    if phi(D) > 0:
        tau1 = mp/pp
    else:
        tau1 = .5 * (1 + mp/pp)
    if phi2(D) > 0:
        tau2 = pm/pp
    else:
        tau2 = .5 * (1 + pm/pp)
       
       
    #if(tau1.imag() < 0):
    #    tau1 = tau1.conjugate()
    #if(tau2.imag() < 0):
    #    tau2 = tau2.conjugate()

    tau1 = reduce_tau(tau1)[0]
    tau2 = reduce_tau(tau2)[0]
    j1 = j_invariant_qexp(10)( CC(2 * pi * i * tau1).exp() ).real();
    j2 = j_invariant_qexp(10)( CC(2 * pi * i * tau2).exp() ).real();

    c4a = (j1 * phi(D)).nth_root(3)
    c4b = (j2 * phi2(D)).nth_root(3)
    if verbose:
        print "Trying c4 = "
        print "  phi(c4) =", c4a
        print " phi2(c4) =", c4b
        sys.stdout.flush()

    # now write c4 as A + B sqrt(5)
    A = ZZ(round(c4a + c4b))/2
    B = ZZ(round((c4a - c4b)/RR(5).sqrt()))/2

    #print "trying c4 around", A, "+", B,"sqrt5"
    base_x = ZZ(round(A - B)) # A - B should be an integer, but
                              # our precision isn't perfect, so it might not be
    base_y = ZZ(2 * B)
    if verbose:
        print "trying c4 around", A, "+", B,"sqrt5 =", base_x + base_y * a
    c4_possibilities = [x + a * y for (x,y) in cartesian_product_iterator([srange(base_x - 5, base_x + 5), srange(base_y - 5, base_y + 5)])]


    #print c4a, c4b
    # check if any nearby polynomials have roots:
    # (it may be better to just compute c4 using c4 = 1/2 trace(c4) + 1/2(trace(c4/sqrt5))sqrt5)
    #r_base = -ZZ(round(c4a + c4b))
    #s_base = ZZ(round(c4a * c4b))
    #c4_possibilities = []
    #for r in range(-3, 3):
    #    for s in range(-3, 3):
    #        print "trying to take c4 as a root of", (y^2 + (r_base + r) * y + (s_base + s))
    #        c4_possibilities.extend( (y^2 + (r_base + r) * y + (s_base + s)).roots())
    return_list = []
    for c4 in c4_possibilities:
        #try:
        #    _ = ZZ(c4)
        #except TypeError:
        #print "Trying c4 =", c4
        Z = c4**3 - 1728 * D
        if Z.is_square():
            c6 = sqrt(Z)
            if verbose:
                print "Finding c6=", c6
            E = EllipticCurve_from_c4c6(c4, c6)
            if verbose:
                print
            E = E.integral_model()
            if verbose:
                print "Found a curve!!!!", c4, c6, E.conductor()
            if E.conductor() == level:
                return_list.append(E)
                #return E
            c6 = -c6
            E = EllipticCurve_from_c4c6(c4, c6)
            if verbose:
                print
            E = E.integral_model()
            if verbose:
                print "Found a curve!!!!", c4, c6, E.conductor()
            if E.conductor() == level:
                return_list.append(E)
                #return E
            if verbose:
                print
        else:
            #print "that c4 is no good"
            pass

    return return_list
